AEDC-TR-66-71 


k 


JUN  H  1966 

APR  7  1967 
OCT  1  i  158E 
APR  1  2  1388 


A  THEORETICAL  REAL-GAS  ANALYSIS 
OF  THE  EXPANSION  TUNNEL 


Glenn  D.  Norfleet  and  F.  C.  Loper 
ARO,  Inc. 


June  1966 


t: eciz^C'IL  ^sports 

FlLz,  CO?** 


Distribution  of  this  document  is  unlimited. 


VON  kArmAn  gas  dynamics  facility 

ARNOLD  ENGINEERING  DEVELOPMENT  CENTER 
AIR  FORCE  SYSTEMS  COMMAND 
ARNOLD  AIR  FORCE  STATION,  TENNESSEE 


PROPERTY  Of  II  o  , 

AS  PORCE 

AF  4 Ot£nr\  1 


References  to  named  commercial  products  in  this  report  are  not  to  be  considered  in  any  sense  as  an 
endorsement  of  the  product  by  the  United  Stales  Air  Force  or  the  Government. 


A  THEORETICAL  REAL-GAS  ANALYSIS 
OF  THE  EXPANSION  TUNNEL 


Glenn  D.  Norfleet  and  F.  C.  Loper 
ARO,  Inc. 


Distribution  of  this  document  is  unlimited. 


AEDC-TR-66-71 


FOREWORD 


The  work  reported  herein  was  sponsored  by  Arnold  Engineering 
Development  Center  (AEDC),  Air  Force  Systems  Command  (AFSC), 
under  Program  Element  62410034,  Project  7778,  Task  777806. 

The  results  of  research  presented  were  obtained  by  ARO,  Inc. 

(a  subsidiary  of  Sverdrup  &  Parcel  and  Associates,  Inc.),  contract 
operator  of  the  AEDC,  AFSC,  Arnold  Air  Force  Station,  Tennessee, 
under  Contract  No.  AF40(600) -  1 200 .  The  research  was  conducted 
from  March  10  to  August  2,  1964,  under  ARO  Project  No.  VJ2447, 
and  the  manuscript  was  submitted  for  publication  on  March  11,  1966. 

The  authors  wish  to  express  their  thanks  to  E.  G.  Burgess,  III, 
C.  H.  Lewis,  and  J.  J.  Lacey  who  were  instrumental  in  the  develop¬ 
ment  of  the  shock  tube  computer  program  used  in  this  investigation. 

This  technical  report  has  been  reviewed  and  is  approved. 

Donald  D.  Carlson 
Colonel,  USAF 
DCS/Plans  and  Technology 


A.  G.  Williams 
Project  Engineer 
Technology  Division 
DCS/Plans  and  Technology 


ii 


AEDC-TR-66-71 


ABSTRACT 


A  theoretical  real- gas  analysis  of  the  expansion  tunnel  is  presented. 
A  digital  computer  program,  developed  for  this  investigation,  is 
discussed,  and  Fortran  listings  and  flow  charts  are  included.  Tunnel 
performance,  test  gas  slug  length,  and  "working"  parameters  are  given 
for  several  expansion  area  ratios.  Driver  temperature  and  energy  re¬ 
quirements  are  given  for  specific  cases. 
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An  aerodynamic  test  device  utilizing  an  unsteady  expansion  to 
accelerate  shock  heated  gas  was  first  proposed  by  Resler  and  Bloxsom 
(Ref.  1)  and  was  treated  briefly  by  He rtzberg  et  al.  (Ref.  2).  Trimpi 
(Ref.  3)  made  a  detailed  theoretical  study  of  the  expansion  tube,  a 
device  in  which  the  entire  expansion  from  the  shock  heated  condition  to 
the  test  condition  is  performed  unsteadily.  In  both  Refs.  2  and  3  it  was 
suggested  that  an  area  change  could  be  added  downstream  of  the  unsteady 
expansion  so  that  part  of  the  expansion  would  be  performed  steadily. 
Trimpi  and  Callis  (Ref.  4)  later  did  a  perfect  gas  analysis  of  such  a  de¬ 
vice,  which  is  called  an  expansion  tunnel.  The  basic  wave  diagrams  for 
both  the  expansion  tube  and  tunnel  are  shown  in  Fig.  1. 

For  a  given  density  level  the  test  gas  velocity  obtainable  with  an 
expansion  tube  is  approximately  twice  that  obtained  by  a  shock  tunnel 
with  the  same  driver  (Ref.  3),  as  shown  in  Fig.  2.  In  spite  of  this 
considerable  performance  gain,  little  effort  has  been  made  to  develop 
an  operable  expansion  tube.  To  a  large  degree  this  hesitancy  comes 
from  the  expected  problems  associated  with  the  device. 

One  of  the  more  questionable  aspects  of  the  expansion  tube  is  the 
uniformity  of  the  test  gas.  Some  insight  into  the  problem  can  be  gained 
from  Fig.  3,  which  gives  the  ratio  of  the  test  gas  slug  length  after  shock 
heating  to  the  tube  diameter.  This  figure  shows  that  the  test  gas  is  in 
close  proximity  to  the  secondary  diaphragm  at  the  time  of  rupture. 
Although  the  secondary  diaphragm  can  be  very  thin,  it  does  require  a 
finite  time  to  rupture.  There  is,  then,  a  shock  wave  reflected  from  the 
diaphragm  which  very  quickly  weakens  to  the  vanishing  point.  Since  the 
diaphragm  will  be  bulged,  this  reflection  will  not  be  planar.  At  least 
the  initial  portion,  and  possibly  all,  of  the  test  gas  passes  through  this 
shock  wave  which  is  not  uniform  in  either  space  or  time. 

For  the  perfect  gas  case  the  addition  of  the  steady  expansion,  as  in 
the  expansion  tunnel,  is  very  effective  in  alleviating  the  test  gas  slug 
length  problem.  *  The  steady  expansion  is  detrimental,  however,  in 
terms  of  test  gas  velocity  obtainable  with  a  given  driver.  2  The  optimum 
expansion  tunnel  design  will  probably  require  some  compromise  between 
maximum  test  gas  slug  length  and  maximum  velocity. 


]See  Fig.  22  of  Ref.  4. 
2See  Fig.  5  of  Ref.  4. 
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The  investigation  reported  herein  was  undertaken  in  order  to  deter¬ 
mine  significant  expansion  tunnel  performance  parameters  based  upon 
a  real  air3  test  gas  and  to  generate  the  proper  working  charts  for 
future  expansion  tunnel  design  and  operation  (Appendix  I). 


SECTION  It 
CALCULATIONS 

2.1  EXPANSION  TUNNEL  CALCULATIONS 

This  investigation  is  concerned  with  the  theoretical  possibilities  of 
an  expansion  tunnel  as  a  high  velocity  flight  duplication4  device.  It  is 
meaningful,  then,  to  determine  the  significant  parameters  such  as  shock 
strengths,  pressures,  and  tube  lengths  in  terms  of  duplication  altitude 
and  flight  velocity.  The  bulk  of  these  computations  was  done  with  a 
digital  computer  program  which  is  briefly  described  below.  A  detailed 
description  including  the  Fortran  listing  and  flow  diagram  is  given  in 
Appendixes  II  and  III.  The  test  gas  slug  length,  it,  and  accelerating 
tube  charge  pressure  were  calculated  by  hand  as  described  in  Sec¬ 
tions  2.  1. 2  and  2.  1.3,  respectively. 

2.1.1  Tki  Expansion  Tunna)  Program 

The  expansion  tunnel  program  is  designed  to  determine  flow  param¬ 
eters  and  tube  lengths  of  interest  for  given  altitude,  velocity,  and 
expansion  area  ratio.  The  input  and  output  data  for  this  program  are 
shown  in  Table  I.  The  expansion  tunnel  portion  of  Fig.  1  shows  the 
various  flow  regions. 

The  general  procedure  is  as  outlined  below. 

1,  For  a  given  flight  altitude,  free-stream  pressure  and  tempera¬ 
ture  are  obtained  by  table  look-up  using  the  data  of  Ref.  5. 
Enthalpy,  hg^,  and  entropy,  Sg^,  are  calculated  using  the 
perfect  gas  equations: 

■  cpTgx 


3 Thermodynamic  and  chemical  equilibrium  was  assumed  for  all 
calculations  included  herein. 

4Dupllcatlon  herein  refers  to  the  complete  matching  of  ambient 
properties  and  ambient  chemistry  together  with  the  required  flow  velocity. 
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'  The  perfect-gas  equations  are  valid  for  altitudes  of  interest. 

2.  The  continuity  and  energy  equations  for  the  steady  isentropic 
expansion  from  6  to  6 A  are  combined,  yielding 

P$A$(2ho6A  *“  2hs)  ;I  =  P§\  1  6A  ^6A 

This  equation  and  the  equation  of  state,  represented  by  the 
Table  of  Thermodynamic  Properties  (Refs.  6  and  7)  are  solved 
(for  S  constant)  simultaneously  for  pg  and  hg. 

3.  The  unsteady  expansion  from  2  to  6  and  the  shock  crossing  1 
to  2  must  be  performed  simultaneously  since  the  limit  of  the 
expansion  is  determined  by  the  shock  crossing. 

For  the  unsteady  isentropic  expansion, 


From  the  tables  of  thermodynamic  properties  (Refs.  5  and  6) 
through  the  unsteady  expansion: 

a  =  (h,S) 


at  2: 


?2  =  h  (1>2  .  s2) 


P2  ~  ^3  (  ^2  »  S2  ) 


Across  the  shock,  Mg  : 


-  1*2  )2  =  Pi  +  Pi  Usj2 

(momentum  equation) 

-  U2)2  =  hx  +  1/2  USl2 

(energy  equation) 

-  u2)  -  Pi  uSl 

(continuity  equation) 

For  the  gas  in  region  1; 

Pi  =  Pi  k]  T i 


(ideal  gas  equation  of  state) 


The  charge  gas  temperature,  Tj,**  and  enthalpy,  hj,  are  inputs 
to  the  program.  The  above  equations  are  solved  for  the  required  flow 
parameters  in  regions  1  and  2,  assuming  the  unsteady  expansion  to  be 
isentropic. 


a 

For  the  calculations  reported  herein,  Tj  =  296°K. 
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The  accelerating  tube  length  is  calculated  from 

18/Atr  =  U6  (M6  -  D 

The  tube  length  ratio  £  is  optimum^  when  the  (U  +  a)  wave 
reflected  from  the  2-3  interface  overtakes  the  tail  of  the  expansion 
at  the  test  section  as  shown  in  Fig.  1.  The  time  required  for  the 
passage  of  this  reflected  wave,  Atw/ig,  is  calculated  by  integration 
through  the  unsteady  expansion.  The  technique  used  to  determine  jCg 
is  described  in  Appendices  II  and  III  and  is  similar  to  that  of  Ref.  8. 

The  theoretical  model  of  this  program  was  based  upon  the  follow¬ 
ing  simplifying  assumptions: 

1.  Air  in  regions  1  and  6A  is  assumed  to  be  ideal. 

2.  Air  in  regions  2  and  6  is  assumed  to  be  in  thermodynamic 
equilibrium. 

3.  Flow  is  inviscid  and  one  dimensional  throughout. 

4.  Diaphragm  rupture  is  instantaneous  with  no  losses. 

5.  The  expansion  nozzle  has  no  length  and  therefore  zero 
"start"  time. 

2.1.2  Test  Gas  Slug  Length 


The  test  gas  slug  length,  was  calculated  from  the  appropriate 
form  of  the  continuity  equation, 

P2  ilA2  =  P6\  USA  Atr  A 6  \ 


2.1.3  Accelerating  Tube  Charge  Pressure 


The  accelerating  tube  charge  pressure,  pg,  was  determined  for  a 
given  ag  using  Ug  and  pg  from: 

Pe  =  P7  (P3/P7)  =  Pe  ( Ps /p7 ) 

and  ,  .  ,  „  . 

P6  =  P7  (across  interface) 

The  ratio  of  Pg/p^  was  obtained  from  Ref.  9  for  a  given 

Ug  /ag  =  U7  /ag  (across  interface). 


6i.e., 


minimum  for  maximum  run  time  with  a  given  jfg. 
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2.2  DRIVER  CALCULATIONS 

The  loss  in  performance  associated  with  the  steady  expansion  can, 
within  limits,  be  offset  by  using  a  higher  performance  driver.  While 
driver  design,  per  se,  is  not  within  the  scope  of  this  investigation,  a 
knowledge  of  driver  requirements  becomes  important  in  assessing  the 
significance  of  the  performance  loss. 

Of  particular  interest  are  driver  pressure,  temperature,  and 
energy  requirements.  Driver  pressure  and  temperature  were  calculated 
using  an  existing  shock  tube  program,  but  driver  energy  (and  driver 
optimization  based  upon  energy  requirements)  was  determined  by  hand 
calculations. 

2.2.1  The  Shock  Tube  Program 

The  shock  tube  program  used  in  this  study  determines  flow  conditions 
for  given  charge  conditions  in  a  two-  or  three-stage  shock  tube.  The 
program  was  used  to  determine  driver  charge  conditions  for  a  given 
and  P|.  This  solution  is  obtained  in  the  usual  manner  of  expanding  the 
driver  gas  to  match  pressure  and  velocity  at  the  2-3  interface.  The 
siiock  crossing  is  performed  much  as  is  done  in  the  expansion  tunnel 
program  except  that  the  thermodynamic  properties  are  obtained  by  an 
empirical  surface  fit  (Ref.  10)  to  the  data  of  Ref.  11.  The  basic  thermo¬ 
dynamic  data  used  in  the  two  machine  programs  differ  in  that  the  data  for 
the  expansion  tunnel  program  include  intermolecular  force  effects  where¬ 
as  the  data  for  the  shock  tube  program  do  not.  The  shock  tube  program 
was  used  only  to  determine  shock  strength,  Mg^,  for  given  driver  and  driven 

tube  conditions.  Experience  has  shown  that  values  of  Mg  calculated  with 
either  s-et  of  data  are  in  good  agreement.  The  inconsistency  in  thermody¬ 
namic  properties  applies  only  to  calculations  involving  the  driver,  and 
there  its  effects  are  felt  to  be  insignificant. 

This  program  will  accept  driver/driven  area  ratios  of  any  value; 
however,  for  <  1  and  low  values  of  p^p  upstream -facing  secondary 
shocks  are  possible.  If  such  a  shock  is  standing  in  the  area  change,  as 
opposed  to  moving  downstream,  the  program  will  not  give  a  solution. 

In  such  cases  the  solution  was  obtained  by  hand  calculations. 

Program  calculations  are  based  upon  the  following  simplifying 
assumptions: 

1.  Gases  are  assumed  to  be  perfect  except  for  the  shock  heated 

region  2  which  is  taken  as  real  air  in  thermodynamic  equilibrium. 
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*-«rr  * 

'•'is 

2.  Flow  is  inviscid  and  one  dimensional  throughout. 

3.  Diaphragm  rupture  is  instantaneous  with  no  losses. 

2.2.2  Driver  Energy  end  Optimization 

Driver  energy,  which  becomes  highly  significant  in  the  case  of  arc 
heating  since  its  magnitude  is  reflected  directly  in  the  cost  of  a  power 
supply,  was  calculated  assuming: 

1.  Helium  as  the  driver  gas  (perfect  gas  assumed). 

2.  The  power  supply  is  of  the  fast  discharge  type,  with  a  discharge 
time  of  the  order  of  100  ^sec.  Arc  efficiency  data  (Fig.  4)  were 
obtained  from  Refs.  12  and  13  and  are  applicable  to  a  fast 
discharge  system. 

Driver  energy  per  unit  volume  for  a  constant  volume  energy  addition 
process  is  given  by 

AE  _  l  P4  T4  -Tj  I 

A4  i4  yA  -  1  i)  T4  **&.• 

Multiplying  by  the  tube  length  and  area  ratios  yields  the  energy  param¬ 
eter,  e, 

=  AE  =  1 _ P4_  T4  -  Tj  /  A6  \  /A4  /  it  ) 

-  A6A  Air  y4  —  l  V  t4  U«a/ \ Aj/WiA  *b/  V  Atr; 


where  pi  and  T4  were  given  by  the  shock  tube  program 

—^4-  was  calculated  by  hand  using  local  values  of  velocity 
1  acoustic  velocity  from  the  shock  tube  program 

— £-e—  and  — l— were  given  by  the  expansion  tunnel  program 

Alr  t  B 


>?  was  taken  from  Fig.  4 
Ti  was  assumed  to  be  296°K 
and  A6  =  A] 


There  are  an  infinite  number  of  combinations  of  T4,  pA,  and  A ^ 
which  will  yield  identical  theoretical  driver  performance.  In  order 
to  give  meaning  to  energy  requirements,  it  is  necessary  to  choose  the 
combination  of  these  variables  which  will  yield  the  given  performance 
while  using  minimum  energy,  i.  e.  ,  the  driver  must  be  optimized.  In 
order  to  optimize  the  driver,  a  power  factor  is  defined  as: 

£  AE  AE  A4  i4 

Aj  / 1  A4  *4  Aj  /  j 
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The  driver  is  optimized  when,  for  a  given  Mg^,  pj,  and  Tj,  the  value 

of  £  (driver  energy  per  driven  tube  unit  volume)  is  minimum.  The 
actual  optimization  process  consisted  of  choosing  values  of  T4  and  A^j, 
calculating  p^,  and  then  varying  T^  and  A^j  until  a  minimum  was 
reached  in  £ . 


SECTION  III 
RESULTS7 


3.1  PERFORMANCE 

Performance  calculations  in  terms  of  velocity  and  altitude  were 
made  for  expansion  area  ratios  of  1,  10,  100,  and  1000.  Required 
shock  strength,  Mg ^ ,  charge  pressure,  p^,  and  compressibility  factor, 

Z2,  are  presented  in  Fig.  5.  The  parameters  Mg^  and  pi  were  chosen 

because  of  their  wide  acceptance  as  independent  variables  in  shock  tube 
work  (particularly  in  shock- crossing  calculations).  Their  values  give 
a  general  indication  of  driver  requirements.  The  value  of  Z2  is  an 
indication  of  the  level  of  dissociation  and  ionization  in  the  shock  heated 
gas  in  region  2.  Bray  (Ref.  14)  has  shown  that  for  the  case  of  a  steady 
expansion  in  which  the  flow  upstream  of  the  expansion  is  not  frozen,  the 
mole  fraction  of  frozen  constituents  after  the  expansion  is  a  very  weak 
function  of  the  ionization-dissociation  level  prior  to  the  expansion.  ® 

If  the  same  phenomenon  occurs  in  an  unsteady  expansion,  then  the  value 
of  Z2  is  not  necessarily  indicative  of  the  ionization- dissociation  level  in 
the  expanded  test  gas.  Until  an  analysis  of  recombination  through  an 
unsteady  expansion  is  available,  no  meaningful  comparison  of  test  gas 
ionization- dissociation  level  for  different  expansion  area  ratios  can  be 
made. 

3.2  DRIVER  REQUIREMENTS 

In  order  to  investigate  driver  requirements  in  more  detail,  driver 
temperature,  T4,  was  calculated  for  a  constant  area  driver  using 
helium  at  a  pressure,  p^,  of  5000  atm.  The  results  are  shown  in  Fig.  6. 
The  drivers  considered  here  represent  rather  severe  requirements,  but 
their  design  is  believed  to  be  within  the  present  "state  of  the  art1'.  The 
high  temperature  helium  drivers  would  probably  be  arc  heated,  although 
the  densities  are  somewhat  higher  than  normal  for  arc  drivers. 


'Additional  results,  in  the  form  of  working  graphs  for  tunnel  design 
and  operation,  are  presented  in  Appendix  I. 

g 

Applies  specifically  to  frozen  atomic  oxygen. 
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The  driver  energy  optimization  was  done  for  a  helium  driver  for 
the  specific  case  of  Mg:  =  10,  px  =  10,  and  Ti  =  296°K.  The  variation 
of  the  power  factor,  £ ,  with  driver  temperature  and  area  ratio  is 
shown  in  Fig.  7.  For  a  maximum  driver  pressure  of  5000  atm,  the 
optimum  occurs  at  a  temperature,  T4>  of  about  5000°K  and  an  area 
ratio  A41  k  1.  The  driver  energy  parameter,  e,  is  shown  for  the 
optimum  driver  in  Fig.  6.  Note  that  these  energies  are  for  constant 
shock  strength  and  charge  pressure,  p^,  and  not  constant  performance. 
Although  they  cannot  be  used  to  compare  energy  requirements  as  a 
function  of  expansion  area  ratio,  they  do  give  single  point  energy  re¬ 
quirements  for  an  optimum  driver  and,  therefore,  some  insight  into 
power  supply  requirements. 


In  order  to  compare  energy  requirements  for  various  expansion 
ratios,  it  is  necessary  to  optimize  the  driver  at  the  same  performance 
level  for  each  expansion  area  ratio.  This  gives  a  different  value  of 
Mgj  and  p-^  for  each  value  of  Ag^/Ag.  In  order  to  simplify  the  optimiza¬ 


tion  it  is  assumed  that  the  optimum  area  ratio,  A^  (0pt}»  eclual  t°  one. 
This  assumption  is  reasonable  for  shock  strengths,  Mg^,  near  10  and 

pressures,  p^,  near  10  atm  since  the  constant  temperature  curves  of 
Fig.  7  are  very  flat  in  the  region  near  A4^  =  1. 


Optimum  driver  temperature  and  pressure  were  calculated  for  the 
performance  level  of  Ug^  =  30,  000  ft/sec  and  altitude  -  150,  000  ft 
(9.  3  <  Mg  <  12.  6,  1. 2  <  Pj  <  10.  4  atm).  Energy  required  for  the  optimum 
case  is  shown  as  a  function  of  the  expansion  area  ratio  in  Fig.  8.  The 
inclusion  of  driver  area  ratio  in  the  optimization  would  produce  second- 
order  changes  in  the  curve;  however,  it  is  doubtful  that  this  would  be 
significant  in  view  of  the  large  variation  of  energy  with  expansion  tunnel 
area  ratio. 


One  point  concerning  the  inviscid  flow  assumption  seems  worthy  of 
mention  here.  There  are  two  viscous  effects  which  can  greatly  affect  the 
driver  energy  requirements: 

1.  Shock  attenuation  (energy  loss) 

2.  Decreased  run  time  (mass  loss). 

In  order  to  offset  shock  attenuation  in  both  the  driver  and  accelerating 
tubes,  a  more  energetic  driver  gas  will  have  to  be  used.  In  addition  a 
longer  driven  tube,  and  therefore  a  longer  driver,  will  be  required  to  re¬ 
cover  the  run  time  lost  by  boundary- layer  effects  (see  Refs.  15  and  16). 

In  terms  of  driver  energy  requirements  the  two  effects  are  additive,  and  it 
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is  likely  that  for  high  SLj d  tubes  the  inviscid  calculation  will  significantly 
underpredict  driver  energy  requirements. 

3.3  RUN  TIME 

Run  time  per  unit  length  of  driven  and  accelerating  tube  is  presented 
in  Fig.  9.  Driver  length  is  not  included  since  it  is  not  defined  by  alti¬ 
tude  alone,  but  depends  upon  the  particular  driver  conditions  chosen. 
Usually,  its  length  will  be  small  compared  to  the  combined  length  of  the 
accelerating  and  driven  tubes.  It  should  be  emphasized  that  the  run  time 
presented  here  assumes  an  instantaneous  nozzle  start.  Although 
Atr/(f  i  +  $. g)  increases  with  increasing  expansion  ratios,  the  effect  of 
expansion  area  on  the  actual  run  time  will  depend  to  a  large  degree  on 
the  nozzle  starting  process.  No  attempt  was  made  in  this  study  to 
determine  nozzle  start  times;  however,  the  ''perfect  start"  perfect  gas 
case  is  treated  in  Ref.  4. 

3.4  TEST  GAS  SLUG  LENGTH 

Figure  10  gives  the  length  of  the  test  gas  slug  at  the  time  of 
secondary  diaphragm  rupture.  The  parameter  on  the  right,  is 

the  ratio  of  the  test  gas  slug  length  to  the  diaphragm  diameter  for  an 
accelerating  tube  length-to-diameter  ratio  of  200.  The  maximum  value 
of  .tfg/dg,  anc*  therefore  is  determined  by  viscous  effects  and  is 

unknown.  However,  from  shock  tube  experience,  anig/dg 
quite  large.  Even  for  the  large  £/d  and  an  expansion  ratio  of  1000,  the 
test  gas  slug  length  is  only  one-tenth  of  the  diaphragm  diameter  for  high 
velocities. 


SECTION  IV 
CONCLUDING  REMARKS 


In  order  to  better  illustrate  the  effect  of  the  steady  expansion,  sum¬ 
mary  plots  (Figs.  11  and  12)  were  made  for  a  specific  case.  For  com¬ 
parison  purposes  the  perfect  gas  results  of  Ref.  4  are  also  shown  in 
Fig.  12. 

Performance  loss  by  the  addition  of  the  area  change  for  an  altitude 
of  150,  000  ft  is  illustrated  in  Fig.  11.  Since  the  performance  lines  are 
for  a  constant  altitude,  and  therefore  a  constant  entropy,  it  follows  that, 
for  a  given  Mg.,  p^  is  constant.  A  given  driver,  then,  would  operate 

along  a  line  of  constant  Mg.. 
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As  shown  in  Fig.  11,  the  loss  in  velocity  for  a  given  driver  can  be 
considerable;  however,  it  becomes  significant  only  when  the  driver  is 
of  limited  potential.  Figure  12a  illustrates  temperature  requirements 
for  a  specific  performance  (Ug^  =  30,  000  ft/sec  at  a  duplicated  alti¬ 
tude  of  150,  000  ft)  and  a  driver  pressure  of  5000  atm.  For  this  per¬ 
formance  point  the  large  expansion  area  ratios  would  entail  severe 
driver  temperature  requirements  for  drivers  which  heat  the  gas  by  heat 
transfer  from  surrounding  surfaces. 

Drivers  which  heat  the  gas  directly,  such  as  electric  arc  heated 
drivers,  have  temperature  limits  which  are  quite  high  and  thus  would 
not  impose  fundamental  limitations  for  the  performance  shown  here. 
Drivers  operating  in  this  mode  are  generally  limited  more  by  energy 
requirements  since  they  add  energy  very  rapidly  and  normally  require 
an  energy  storage  system  which  has  a  high  relative  cost.  For  a  given 
test  section  size  and  run  time,  driver  energy  and,  therefore,  power 
supply  costs  decrease  with  increasing  expansion  area  ratio  (Fig.  12b). 

For  a  given  accelerating  tube  length-to-diameter  ratio,  J?g/dg, 
large  gains  can  be  made  in  test  slug  length,  f^/dj  (and  hopefully  flow 
uniformity),  by  using  large  expansion  area  ratios  (Fig.  12c).  It 
should  be  noted  that  the  parameter  actually  of  interest  is  it/d^  (max) 
which  occurs  at  ig/dg  =  ig/dg  (max)’  anc*  dependent  upon  boundary- 

layer  growth.  No  attempt  is  made  here  to  include  boundary-layer 
effects;  however,  perfect  gas  calculations  for  a  simplified  model  are 
included  in  Ref.  4. 

In  summary,  an  increasing  expansion  area  ratio  causes: 

1.  A  loss  in  performance  if  the  driver  is  limited  in  temperature 
and  pressure. 

2.  A  gain  in  performance  if  the  driver  is  limited  only  in  energy. 

3.  An  increase  in  the  test  gas  slug  length  parameter,  At/#g- 
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Fig.  1  Wave  Diagrams  -  Expansion  Tube  and  Expansion  Tunne 
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Fig.  2  Comparison  of  Expansion  Tube  and  Shock  Tunnel 
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Fig.  4  Arc  Driver  Efficiency 
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o.  A^a/A^  =  1 
Fig-  5  Performance  Map 
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Fig*  5  Continued 
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Fig,  6  Driver  Requirements 
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c.  &6A  “  1  00 

Fig.  6  Continued 
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Fig.  10  Proximity  of  Test  Gos  to  Secondary  Diaphragm  —  Expansion  Tunnel 
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c*  Test  Gos  Slug  Length 

ig.  12  Variation  of  Some  Critical  Parameters  with  Expansion  Area  Ratio 
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TABLE  I 

EXPANSION  TUNNEL  PROGRAM  -  INPUT  AND  OUTPUT  DATA 

INPUT  DATA 

U6A  Ft  (Altitude)  A6A/A6  T1 


OUTPUT  DATA 


pl 

us 

bl 

S1 

^l/ig 

V“r 

P2 

U2 

T2 

P  2 

h2 

a2 

Z2 

p6 

U6 

T6 

p6 

h6 

a6 

Z6 

P6A 

U6A 

T6A 

p6A 

h6A 

a6A 

Z6A 

30 


A  EDC-TR-66-71 


APPENDIX  I 
WORKING  GRAPHS 


As  an  aid  in  expansion  tunnel  design  and  operation,  some  of  the 
more  meaningful  parameters  are  presented  in  the  form  of  working 
graphs.  Charge  pressures  and  pg  are  presented  in  Figs.  1-1  and 
1-2,  respectively.  Pressure  in  the  shock  heated  region,  P2,  is  pre¬ 
sented  in  Fig.  1-3.  The  nondimens ional  form,  p/pg^,  reduces  the 
variation  with  altitude  to  that  caused  by  real -gas  effects  and  acoustic 
velocity  variation. 

Shock  strengths  as  a  function  of  altitude  and  test  gas  velocity  are 
presented  in  Figs.  1-4  and  1-5.  Here  again,  the  variation  with  altitude 
is  caused  by  variation  of  acoustic  velocity,  ag^,  and  real- gas  effects. 

The  optimum  driven  tube  length  is  given  as  a  function  of  test  gas 
velocity  and  altitude  in  Fig.  1-6.  As  noted  previously,  viscous  effects 
in  the  driven  tube  may  increase  the  optimum  length  considerably. 
Accelerating  tube  length  per  unit  run  time  is  shown  in  Fig.  1-7.  The 
effects  of  viscosity  on  run  time  have  not  been  assessed,  even  qualita¬ 
tively,  to  date. 
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x  10'^  ft/sec 


a.  A^a  /A&  1  *0 

Fig,  1-1  Charge  Pressure  —  P] 
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b.  Afix/Afi  =  10 
Fig.  1*1  Continued 
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c •  A$a  /A  g  -  100 

Fig,  [-1  Continued 
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U6A  x  10-3  ft/sec 


d.  A$a  6  =  1  000 

Fig.  1-1  Concluded 
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U6A  x  10"^  ft/sec 

b.  A6A/A6  *  1000  ond  A6a/A5  «  100 
Fig.  1*2  Concluded 
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a.  A 6A  /A 6  =  1.0 

Fig.  1-3  Pressure  behind  Primary  Shock  in  Driven  Tube  —  P2 
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x  1CT3  ft /sec 


c.  A$a  /Afi  “  ICO 
Fig.  1-3  Continued 
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d.  A6a/A6  -  1000 

Fig.  1-4  Concluded 


Altitude  ■  100, 000  -  250, 000  ft 
Air  in  Region  8 
T8  =  296°K 
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APPENDIX  II 

DERIVATION  OF  EQUATIONS  USED  IN  THE  COMPUTER  PROGRAM 
SHOCK  CROSSING  EQUATIONS 

For  a  plane  shock  wave  moving  into  a  quiescent  gas,  the  equations 
for  the  conservation  of  mass,  momentum,  and  energy,  and  the  equa¬ 
tion  of  state  are,  respectively. 


P2  («S,  ~  U2)  =  Pl  USj  (II-  1 ) 

P2  +  P2  (Us,  -  U2)2  =  Pl  +  Pj  USl2  {II- 2) 

h2  +  1/2  (USl  -  l2)2  =  hi  +  1/2  USl2  (II-3) 

P,  =  P,  RT,  (II- 4) 


These  four  equations  will  be  used  to  determine  pp  Ug^,  p^,  and 

2 

Eliminating  U2  and  Ug^  from  Eqs,  (II- 1),  {II- 2),  and  (II- 3)  gives^ 

(P2/Pl  -  l)  -  (11-5) 

2  1  Pl  1  +  Pl/P2 


Now,  eliminating  pj  between  Eqs.  (II-4)  and  (II-5),  one  gets,  after 
some  manipulation. 


P 1  +  P 


2P2  (^2  ~  ^1 )  _  P2  T«  "I  _  P2P2 

1  |_  RT  j  P2  T,  J  Tx 


=  0 


(II- 6) 


The  positive  sign  in  the  quadratic  formula  corresponding  to  the 
above  equation  yields  the  desired  value  of  pj.  Equations  (II- 3),  (II-4), 
and  (II- 1)  can  be  written 


u*i 


(II- 7 ) 


“Introducing  the  nondimensional  quantities  of  p  in  atm  and 
p  in  amagats. 
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Pi  =  Pi  Ti  / Ta 


(11*8) 


U2  =  USl  U  -  Pi/P2]  (n-9) 

Use  is  made  of  Eqs.  (II- 6),  (II- 7),  (II- 8),  and  (II- 9 )  in  Appendix  III. 

LENGTH  EQUATIONS 

Equation  (II- 3),  also  appearing  in  Appendix  III,  will  now  be  de¬ 
rived.  Reference  to  Fig.  II- 1  will  aid  in  the  understanding  of  the 
following  development.  Hence,  from  Fig.  II- 1  it  follows  that 


ii/Aty  =  a2  USl/(USl  -  u2> 
ift/Atj.  =  U6  -  a6 


(II-  10) 
(II- 11) 


For  the  unsteady  expansion  between  regions  2  and  6,  the  following 
equations  are  valid: 


dt 
d  t 

l  -  Al* 


1 


U  +  a 
i 


t  U  -  a 

adU  =  -dh 


(11-12) 

(11-13) 

(11-14) 


The  parameter  i  can  be  eliminated  between  Eqs.  (11-12)  and 
(11-13)  by  differentiating  (Eq.  11-13)  and  putting  the  result  into  Eq.  (11-12). 


(U  -  a)dt  +  (dTJ  -  da )  ( t  -  Atx)  =  (U  +  a)dt 


(11-15) 


Using  Eq.  (11-14)  to  remove  dU  from  Eq.  (11-15)  and  simplifying 


gives 


2dt  da 


t  -  Au 


dh 


(11-16) 


Integrating  Eq.  (11-16)  between  regions  2  and  6  gives 


Hence, 


h2 

log 

a?  ( ~  Atx 

=  -  /  -^-  +  log 

a6  (t6  -  At*  )2 

^6  a 

sum 

a. 

(11-17) 


52 


AEDC-TR-66-71 


where 


so  that 


sum 


aa 

r  a  T*a 


J 

b6/R 


b  2/  R 


d  (  h/  R  ) 
(a/a.  )2 


Aty 

"At” 


exp  (-  0.5  sum/aa) 


(11-18) 


Combining  Eqs.  (11-10),  <11- 11),  and  (11-18)  gives 
<i  =  _ i 

*  0  ,  Ue.  -  U2  /ao  \  '2 

( K 6  *■  a g  J  i - ( — )  exp  (0.5  suan/aa  ) 

8 2  US,  '  a6/ 

or,  in  the  terminology  of  Appendix  III, 


where 


o  _  1 _ 

1  (  U  g  —  a  g  )  y 


>'  =  — - — — (")  exp  (0.5  sum/aa) 
®2  l-S,  V  a6  / 


i8  =  1 
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The  symbols  associated  with  the  above  lines  represent 
the  reciprocal  of  the  slope  of  the  lines. 


Fig.  11-1  Wave  Diagram  Illustrating  Nomenclature  Used  in  Calculation  of  Tube  Length 
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APPENDIX  III 
COMPUTER  PROGRAM 


The  purpose  of  this  program  is  to  compute  certain  quantities  of 
interest  associated  with  the  unsteady  expansion  problem  among  which 
are  the  6  conditions,  the  2  conditions,  and  a  characteristic  length. 

The  method  by  which  the  above  is  accomplished  basically  reduces 
to  solving  two  separate  nonlinear  equations  by  iterative  schemes,  and 
evaluating  certain  integrals  by  Gaussian  mechanical  quadrature.  The 
gas  properties  are  retrieved  from  magnetic  tape  using  a  double  four- 
point  interpolation. 

Program  flow  is  as  follows  (Figs.  Ill- 1  and  2,  and  Tables  III-  1 , 

-2,  -3,  and  -4): 

A.  SR,  if  not  read  in,  is  computed  from  T6A  and  P6A,  TGA  and 
PGA  being  computed  as  functions  of  altitude  in  subroutine  ATP. 

B.  The  remaining  6A  conditions  are  computed  from  SR  and  P6A. 

C.  An  iteration  is  performed  to  find  p6.  The  initial  guess  is 

p  6  =  p  64  -MA. 

A6 

Given  a  value  of  p6,  H6  is  computed  as  a  function  of  SR  and 
p6  (table  look-up).  The  value  of  p6NEW  is  defined  as 

p6NF.lt  - 

v  2<H06A-H6>  A6 

If  p6  and  p6NEW  differ  by  not  more  than  0.  01  percent,  then  the 
iteration  is  said  to  have  converged.  Otherwise,  the  iterative 
value  of  p6  is  taken  to  be  the  average  of  p6  and  p6NEW,  and 
the  process  is  continued. 

D.  The  remaining  6  conditions  are  computed  from  SR  and  H6  by 
interpolation  in  the  tables. 

E.  An  iteration  is  performed  to  find  H2.  The  initial  guess  on  H2 
is  taken  to  be  10 

H2  =  0.6H1  (1  +  0.195X2) 


^Developed  empirically  from  the  perfect  gas  results  of  Ref.  4, 
page  60. 
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where 

X  =  0.32  -IffiA  +  1 

A6A  y  A 6  J 

For  known  H2,  U2  can  be  computed  by  two  relations.  The 
object  of  this  iteration  is  to  find  the  value  of  H  such  that 
these  two  values  of  U,  U21  and  U22,  are  equal.  The  two 
relations  that  U  must  satisfy  are 


P  2  - 

■  p 2  (SR, 

H2  )  (from 

tape) 

(iii- 

la) 

P2  = 

=  P2  (SR, 

H2  )  (from 

tape) 

(iii- 

lb) 

B  = 

2  p  2  (H2  -  HI)  t 
RT1 

P2  -  : 

"(#) 

(m- 

lc) 

C  ■ 

*  -P2p2 

TA/Tl 

(iil- 

Id) 

pi  - 

=  0.5  (-B 

+  s/W 

-  4C) 

J  / 

(iii- 

le) 

US1  • 

-  ^2(H2 

-  m  )/(i 

Pi* 

pi2 

)) 

(iii- 

If) 

PI  . 

=  pi  (Tl/TA) 

(III- 

lg) 

U21  . 

=  US1  a 

-  pl/p2) 

(in- 

lh) 

H2/R 

U22 

-  U6  - 

are  f 

/ 

d (H/R  ) 

(III- 2) 

rA  TA  H 

67  R 

a/are  f 

liquations  III- 1  e  through  III- lh  correspond  to  Eqs.  11-6  through 
II-9,  respectively.  Notice  that  all  the  quantities  necessary  to 
compute  U21  and  U22  are  known  for  an  assumed  H2.  Newton's 
method  with  numerical  first  derivative  is  used  to  find  the  root  of 

f(H2)  =  U22  -  U21  =  0 

Hence 

where 


As  this  technique  requires  two  initial  guesses,  the  second  initial 
guess  is  taken  to  be  1.1  times  the  first  initial  guess.  With  the 


H2i+1  =  H2,  - 


HH2,) 
g(H2, ) 


g(H2, )  - 


Jf  CH2|  )  -  f  (  H  2;  —  i  ) 
H  2,  —  H  2,  —  ! 
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indicated  initial  guesses  and  iteration  method,  the  root  has 
converged  in  every  case  thus  far.  The  iteration  is  said  to 
have  converged  whenever  H2^  +  j  and  H2^  differ  by  not  more 
than  0.  01  percent. 

F.  The  remaining  region  2  conditions  are  computed  in  terms  of 
SR  and  H2. 

G.  The  characteristic  length,  LI,  is  obtained  as  follows: 


rsi  -  ua 

a2  US1 


<III-3a> 


SUM  =  / 

d ( H/R  ) 

(III- 3b) 

r  T  » 

*  La  H6/R 

{  a /a  ref  )2 

Y  ■  “  (S)* 

exp  (0.5  SUM/are £ ) 

(III-3c) 

11  -  1 

(U6-a6)Y 

(III- 3d) 

Input  to  the  program  is  read  from  two  different  tapes. 

A.  Tape  J1N2  (J1N2  =  10) 

This  tape  contains  the  gas  properties  mentioned  above. 
See  Subroutine  SLOW  for  the  proper  format. 

B.  Tape  J1N1  (J1N1  =  5) 

Two  read  statements  are  executed  by  this  tape. 

1.  A  title  card.  Format  (72H.  .  .  ) 

2.  Input  data.  Format  (6E12.  0,  12) 

a.  U6A  (ft/sec) 

b.  T1  (TO 

c.  A6AA 

d.  SR 

e.  P6A  (atm) 

f.  ALT  (ft) 

g.  MORED 

All  output  from  the  program  is  on  tape  JOTJT  (JOUT  =  6). 

A.  The  title  card 

B.  SR 
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C.  Inputs  and  other  constants 

1.  Tl 

2.  HI 

3.  A6A/A6 

4.  RHOA 

5.  PA 

6.  R 

7.  SPEED  REF 

D.  Values  of  the  following  at  the  6-A,  6,  and  2  conditions. 

1.  P 

2.  U 

3.  T 

4.  p 

5.  H 

6.  A 

7.  Z 

E.  Other  output 

1.  LI 

2.  PI 

3.  US1 

4.  MSI 

5.  DT 

6.  L8/DTR  (for  unit  L8) 


Subroutine  INTRP 

The  purpose  of  this  routine  is  to  do. an  N-point  Lagrange  inter¬ 
polation  where  N- 1  is  a  natural  number.  The  argument  list  is: 

(N,  X,  Y,  XINT,  YINT) 

N  is  the  number  of  points 

X  is  the  set  of  independent  values  ^ 

Y  is  the  set  of  dependent  values^ 

XINT  is  the  value  of  the  independent  variable  at  which 
the  interpolation  is  to  take  place 

YINT  is  the  interpolated  value  of  the  dependent  variable 
(the  return  argument) 


11 


X  and  Y  should  be  appropriately  dimensioned  in  the  calling 


routine. 
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Subroutine  GAUSS 


This  subroutine  defines  constants  and  x^  (i  -  1,  16)  such  that 


/  f(x)<k 

? 


can  be  approximated  by 


q  -  p 


16 

V 


bif 


(Xi 


q  r  p 


The  values  of  b^  and  were  taken  from  Ref.  17. 

The  argument  list  for  this  subroutine  is  (b,  x).  Both  b  and  x  should 
be  dimensioned  sixteen  in  the  calling  program. 

Subroutine  SLOW 

The  purpose  of  this  subroutine  is  to  do  a  cross  four-point  central 
Lagrange  interpolation  of  data  which  have  been  stored  on  tape  as  a 
function  of  two  independent  variables.  The  manner  in  which  the  input 
tape  has  been  created  should  be  equivalent  to  the  following: 

DO 1  1  K  =  l,  N 

1  WRITE  (IT)  X  (K),  J,  ( (Y  ( K,  I,  L),  I  =  1,  NV ),  L  -  1,  J) 
where  4  <  N 

2  £  NV  <  9,  a  constant  defining  the  number  of  variables 
exclusive  of  X  and  J 

4  £  J  £  150,  a  variable  defining  the  number  of  points  for  a 
given  K 

For  a  given  value  of  K,  it  is  required  that  Y(K,  I,  L)  be  a  strictly  mono- 
tomic  function  of  L  for  at  least  one  I.  It  is  also  required  that  X(K)  is  a 
strictly  monotonic  (increasing  or  decreasing)  function  of  K. 

While  X  must  always  be  one  of  the  independent  variables,  the  second 
independent  variable  and  the  dependent  variable  need  not  be  specified 
until  call  time.  Any  Y  that  is  a  strictly  monotonic  function  of  K  can  be 
used  as  the  second  independent  variable. 
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The  meaning  of  the  variables  on  the  tape  associated  with  this  par¬ 
ticular  problem  is  as  follows: 


Fortran  Name 

Identification 

X(K) 

SR 

J 

variable 

NV 

9 

Y(K,  1,  L) 

T 

Y(K,2,  L) 

log,0  (p/pA ) 

Y (  K,  3,  L) 

log10  (P/PA) 

Y  ( K,  4,  L) 

log10  (H/R) 

Y (K,  5,  L) 

ye 

Y(K,  6,  L) 

a/aref 

Y (K,  7,  L) 

Z 

Y  ( K,  8,  L) 

H/RT 

Y (K,  9,  L) 

Z* 

The  data  on  this  tape  were  taken  from  Refs.  6  and  7  primarily;  however, 
certain  unpublished  extrapolations  of  the  above  are  also  present. 

The  argument  list  for  the  subroutine  is 

(XX,  Z,  II,  Jl,  IT,  NV,  NERR) 

XX  is  a  specified  value  of  SR 

Z  is  a  subscripted  variable  dimensioned  appropriately 
in  the  calling  routine 

II  is  a  subscript  indicating  that  Z{I1)  is  the  second 
independent  variable 

J]  is  a  subscript  indicating  that  Z(J1)  is  the 
dependent  variable 

IT  indicates  the  channel  and  unit  number  on  which  the 
input  gas  tape  is  mounted 

NV  indicates  the  number  of  variables  on  the  tape  corres¬ 
ponding  to  a  value  of  XX  (NV  is  nine  in  this  case) 

NERR  will  be  returned  equal  to  one  if  and  only  if  the 
interpolation  failed  for  any  reason. 

Subroutine  ATP 

This  subroutine  was  not  written  specifically  for  this  program,  and 
hence  has  options  not  used  here.  Use  is  made  of  the  subroutine  in  this 
program  to  find  temperature  and  pressure  as  a  function  of  altitude.  The 
data  used  by  ATP  were  taken  from  Ref.  5. 
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Subroutine 

GAUSS 


Main 

Program 


l _ , 

_ i 

Subroutine 

ATP 

Fig.  Ill-l  Tree  Diagram 


_ i 

Subroi 

SL( 

. 

utine 

)W 

Subroutine 

INTRP 
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Fig.  111-2  Flow  Chart  of  Main  Program 
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TABLE  HI-1 

CROSS  INDEX  OF  NOMENCLATURE 

As  Used  in  Appendix  III  As  Listed  in  Report  Nomenclature 


ALT 

SR 

S/R 

HO 

ho 

H06 

ho6 

TA 

Ta 

yA 

Ya 

are! 

Qa 

R 

R 

HI 

hi 

Ti 

Ti 

Pi 

Pi 

PI 

Pi 

US1 

Usi 

H2 

h2 

T'2 

t2 

p2 

P  2 

P2 

P  2 

U2 

U2 

H6 

h6 

T6 

T6 

p6 

P6 

P6 

Pe 

U6 

u6 

T6A 

t6. 

p6A 

P6. 

U6A 

U6a 

P6A 

P6a 

A6A 

a<5. 

a 

a 
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TABLE  III •  2 
SAMPLE  INPUT 


Date  Set  1 

(A)  Card  1  (Title  card) 

(B)  Card  2 

1.  U6A  =  20000. 

2.  T1  =  296. 

3.  A6AA  =  10. 

4.  SR  =  (Not  required) 

5.  P6A  =  (Not  required) 

6.  ALT  =  150,000. 

7.  MORGD  -  0  (indicates  another  set  of  data  follows) 


Data  Set  2 


(A) 

Card 

1  (Title 

(B) 

Card 

2 

I. 

1J6A 

2, 

T1 

3. 

A6AA 

4. 

SR 

5. 

P6A 

6. 

ALT 

7. 

MORLD 

ard) 

=  8020. 

=  296. 

=  1. 

=  26.62 
=  .9 

=  (Not  required) 

-  1  (indicates  end  of  job) 
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TABLE  111-3 
SAMPLE  OUTPUT 

SAMPLE  INPUT  Shot  one  UN51E40Y  EXPANSION  IUPUT  ALT ITUDE 1 
5/R-C. 3Q10a0B7E  C2 


I  NPUT 


n 

0.2960COO0E  D3 

HI 

0.319979m 

07 

A6A/A6 

0.099999996 

02 

RHOA 

0.099999996  01 

PA 

0.09999999F 

01 

R 

0. 309170006 

04 

SPEED  REF 
0.10873999E  04 

6-A  CONDITIONS 

P 

0.13428957E-C 2 

u 

0.200000006 

C5 

T 

0.2&564383E 

03 

RHO 

0.138027706-02 

H 

0 » 2865720  3E 

07 

A 

0.107280736 

04 

2 

0,999999946  00 

6-COND I  T  1 0 Kl S 

P 

O.33525464f-Ol 

U 

0.197848776 

05 

T 

0. 635980626 

03 

BHO 

0.139528486-01 

H 

0.71796305E 

07 

A 

0.166&B6&1E 

04 

l 

0.99999993E  00 

2-CONOIUONS 

P 

G . 2  2  4  L4427E  02 

U 

0*735791786 

04 

X 

0,293141906 

tH 

PHD 

0.207423846  01 

M 

0.392045646 

o& 

A 

0 .334448686 

04 

l 

0.10069432E  01 

Ll.. 750296106-01 


PJ-. 32025371E-00 
USl«.B5ftG4500£  04 
NSL-.758U1096E  01 
0T- »  735626346-03 
l fl/OTH-.21S04901E  06 


SAMPLE  INPUT  SHQ7  TNC  UNSTEAQV  EXPANSION  I  INPUT  5/R.P/Pa) 
S/ft-O. 26620000E  02 


INPUT 


T 1 

0.29600000E 

03 

HI 

0.3 199T999E 

07 

A6A/A6 

0.099S9999E 

01 

RHOA 

0.09999999E  01 

PA 

0.09999999E 

01 

R 

0. J0917000E 

04 

S^EFO  REF 
0.108T3999E 

04 

6-A  CONDITIONS 

9 

C.90000000E  00 

It 

0.BO20QQ0OE 

04 

7 

0.62182069E 

03 

RHO 

0*394934046— CO 

H 

C. 6790091 16 

07 

A 

0. 1626098S6 

04 

z 

0 ♦ 10003  L  26E 

Oi 

6-C0ND1TICNS 

P 

0.900066646 

OO 

U 

0.80 199 99 9E 

04 

T 

0,621843376 

03 

PhO 

0.39493404E-00 

H 

0.679015706 

07 

A 

0.16261 16  3E 

04 

Z 

0.10003123E 

01 

2-C0NO1T IQN5 

P 

0. L4919058E 

02 

U 

0 . 39505739E 

04 

T 

0.127883  L  5E 

04 

RHO 

0.317284046  01 

H 

0. 14T95274E 

03 

A 

0.23004321E 

04 

l 

0.10043706E 

01 

L1..5443B253E  00 


P l ■ .672 1 0025E  CO 
US1«=.49104284E  04 
NSI-.43379S25E  01 
DT-.ZB719155E-03 
L0/OTR*. 31534 56 2E  05 
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TABLE  111-4 
FORTRAN  LISTING 


*]p  LOPEP  IBJ08  UNSTEADY  EXPANSION 

$JOBOPOMAPiLOGIC 

siipap  uNrrs 

FNTRY  -UN  10. 

.U**]0.  PI E  UNIT10 

UN  I  T 1 0  file  *Bl  1  1 -iREADY.  INOUT  .8LK32S6.L0W.BIN 

END 

J1BFTC  MJkENP  M94.XR7 

C  FORTRAN  4  UNSTEADY  EXPANSION  PROGRAM 

C  SUBROUTINES  REQUIRED  ARE  GAU5S. ATP. SLCvN  [NTRP 

DIMENSION  U  7  l  »B(161  »wn6J 
COMMON  IMFT 
CALL  CAUSS  I  8  » V 1 
RHOA=  1* 

PA*  1. 

PAA*  14.695 
R*  2001.7 
ASPEO  =  10B7.4 
GAMA=  l.A 
TA=  273.15 
CP1=  10810.13* 

C 

c  tape  order  FOR  CONSTANT  S/fi 

Cl  J  .NUMBER  V.F  POINTS  FOR  THIS  S/R 

C2  r  .TEMPERATURE 

C3  alog IPHO/RHOA | 

C  4  ALOGIP/PAJ 

C5  ALOGIH/Rl 

C6  GAME  i GAMMA  5U3  E 

C  7  A/AA  A/CA  SUB  A  ) 

ce  z 

r«?  h/rt 

CTO  Z*  l  STAR 

c 

5  F0RKATf72^ 

1  > 

16  FORMAT i  6H0  S/R*€ 14 « 8 / / ?H?  INPUT/114H  T1 

1  A6A/A6  RhOA  PA  R 

2  SPEED  REF  /1H  7F16.8  6-A  CONDITIONS  /114H 

3  P  U  T  RhO 

4  H  A  2  /1H  7E16-8/1 


HI 


2S  FORMAT  C  ]  *H0  6-CONDlTlONS  /11<iH  P  U 

3  T  RHO  H  A 

2  Z  /1H  7E16.8/  1 

64  FORMAT  USHO  ?-CONDlTlONS  /114H  P  U 

1  T  RHQ  H  A 

?  Z  /1H  7E16.8/  } 

91  FORMAT <5H0  Ll*E13.e/5H0  Pl*CI3  . 8/6H0  USl^En.e/fcHO  MS  1 =E 1 3 . 6 /5hO  0 
lT=E12.6/9Hn  L0/OTRcE12.G  ) 

C  J ] N 1 =  INPJT  OAT  A  TAPE 

J1N1=S 

C  J | N 2 s  INPUT  GAS  PROPERTY  TAPE 

JIN?=10 

C  JOUT*  OUTPUT  DAI  A  TAPE  I  THE  ONLY  CUE) 

JOUT=b 
NT1MF»1 
1  C0NT1NUF 
TMFT.O 

4  READ  IJ1NI.5) 

WRITE  I J0UT*5> 

REAP  CU1N1.2J  U6A*TI.A6AA.SR.P6ApALT.M0RED 
?  FORMAT  t  6  F  3 ?.C*1?  > 

Hi  =  CP  1  *  T 1 

1CC ALT) 200*10. 2C0 

200  CALL  ATP  C 1 *2  * t .2.ALT . T6A.P6A.DMH .DUMJ 
P*A=P6A/PAA 

SR=  23. =<56  +  2 . 5+ A  LOG  1  T  6  A  /  T  A }  -ALOG  I  P6A/P  A  I 
]0  NFftR-D 

2  C  3  =  AL0010  CP6A/PA  1 

CALL  SLOW  f SR.Z»3»1 .J1N2.9.NERR  ) 

[FINERR-im  .IOC.]  1 
T 1  T6  ft  =  2  t  l  ) 

CALL  SLOW  1 Sft,2 .3 .2 rJINJ.O.NERR  I 
12  PH06A-  RHOA^l G.««Z C23 

CAuL  SLOW  fSR.Z.B.A.JIN.’^.fiERR  > 

M6A  *  g«]0-“«ZUl 

CALL  SLOV  <SR.2,3,6 i JIN?r9.NERR  ) 

14  AftA=  A&PE0«£I6I 

CALL  SLOW  ISR.2 .3.7.J1N2.9.NFRR  ) 

|S  2* A c  H  7) 

m064=  MAA*.5*U6A  "LriA 

WRlTF  ( JOUT i  1 6  I  SftiTI  .Hi  iA6AA»RHOA.P4.R,ASPEn.P6A.UtA.T6A. 

1 PH06A,H6AiA6A .ZftA 
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TABLE  1 11-4  (Continued) 


RHC6  =  RH06A* A6  AA 

DO  18  r»l**OQ 

112)=  4LOG1DI RhO*/RHCA» 

CALI,  SLOW  <  SR*Z*2»4*JtN2*9*NERR  I 
tF(NFRR-n  17* ICC  <  1  7 

IT  H6  = 

TE'-IP-  5CPT4 2**(HQ6A-H6) ) 

TFMPx  RH06A*U6A*A6AA/TEMP 
IFC  A*IS  I  tT^MP-RHOfeJ /TEMPI  MOO. -.01  >20*20*]  ft 
18  RhO&=.5« (RH06+TEMP > 

1?  PAUSE  nil] 

GO  T0  IOC 

20  CALL  SLOW  {SR*Z*4*1*J1S2*?*NER8  3 
U6  =  RH0FA*UfcA*A6AA/RHC6 

If 1NERR“11?1*100»21 

21  T6*2<U 

CALL  SLOW  |SR*Z*4*3*JIN2*9*NERR  > 

22  PS  =  PA*10**+Z < 3 ) 

CALL  SLOW  [SR*Z*4.6*J[N?.9,\ERR  I 
28  AA=  ASP'^ZIS) 

CALL  SLOW  <5R»2i*'YtJ(N?«9f*ERR  > 

24  Z6=Z(7] 

WR JTE  t JCUT*2S>  P 6  * U6 .  T6 • RK06  *H6  *  A6  » Z  6 
GO  TO  30 

C  subroutine  TO  COMPUTE  INTEGRAL 

57  5UM=0, 

1  =  0 

2e  r=i+i 

2<a)  =  .MCrfl  !  [*(  XU-XL  l+XU  +  XLl 
l\ 4  I  *  AL0G10 \ Z  I u  u 

CALL  SLOW  ($R*Z,4i6-*JtN2.9.NERR  ) 

I p  (NERR-1 129. 100.29 
2  9  SUM.  SLIM  ( I)  / <  Zl6l**KC0  > 
tFU  16128*26*26 
26  COaTINLF 

SUM.  SC V *  I  XU-X4.1**S*AS^E0/<GAMA*TA> 

GO  TO  (3A.76l.KGO 
r  the  ^-conditions 

?0  XX  =  • 32*U6 A/ A6A  + ( A6AA > ** . 26  + 1  - 

M2>*  H'*M  M.  19S*%X»XX  l«.6 
H2=H2 1 


00  SO  J=1,S00 
Z<4|-  ALOG 1 0 ( H2 / R  3 
CALL  5 LOW  (SFs.Z  ♦4*2»U1N2.9*NFRR  [ 
iFiNFRR-mi.ioo.il 
81  RH02=  RHOA* 1 0 • * *Z ! 2  ) 

call  Slow  [Sr*z*4***jim2*«iner<?  ) 

^2  P2 =  P  A* 1 0*  **2 I3| 

PB=  (2*»»RH02*<H7’HI  |/<R*T1  |J*RHC2  P2*TA/T1 

CC=  -P2*RH0?*TA/T 1 

RH01  =  *S+I -3B+SCRT  rfiBMB-A.KC  J  J 

USl  =  SCRT1  I  f<  1  .  Rh01«PH0. /<  RHC2MH0?)  ») 

P!  sRMOl *T  1  / T  A 

L22  =  USl  *  I  1  .-P.H01  /RH02  I 

XU=  H2/=? 

XL=  HG/R 
KGC=  ] 

GO  TO  27 
?5  U2  !  =  UA-Sr  M 

!r  t*>-\  f4i> 

'6  FX1 =  U21-U2? 

H?2  = 

GO  TO  A R 
AO  FX2  =  U21-U22 

RAT  JO-  I  *-H21/H221 / I 1.-PX1/FX2) 

H?3-  H22* I  1 *-7AT 10 1 

I F  (  Ac  5  I  ( -I? A-H27  )/H2 3)  MOD. --0J  )6D*60*a«5 
t-S  FX1-FX2 
H21 *H?2 
H22=H2  * 
a 9  H2 =  H2? 

*0  CCNT[N'JF 

PAU5F  ??227 
SC  TO  ICC 

6C  2C43>  AL0G10CH?t/R) 

CALL  SLOW  ISR.2 ,A*1*JTM2.9.NERR  i 
|F<NFRR-1 )61 *100.61 
61  T?*Zl 1  I 

CALL  SLOW  f SR. 2 . A *6* JIN2*9 »M~R3  ) 

6?  A?=-ASCED»7»6' 

CALL  SLOW  :SR*Z  .4  i  7  ♦  J  1  N2  •<?  .  NE  PR  I 
63  22=2(7) 

U?=U?1 
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TABLE  111-4  (Continued) 


VR] 7* < JCUT.641  P?  U2*T2,RH0?.H2*A2.Z? 

C  characters 1C  length 

7C  OFLT.  A2*U5)1 

Y*  DFLT 
*U  =  H?/F> 
xi~  H6/R 
KGQ-? 

GO  TO  27 

7  6  Y=Y*S09T< a?/A6 J*E>PI . 5* SJM/ASPFD) 

XL  1  =  1  * /  I  <U6-A6>»Y J 

X*S1  =  USJ*S0RT ITA/T]  > /ASPEO 

XLf.nT=  U<>*<U6-Afc>/Ai 

W^lTF  I  JOiJ T  *9  1  1  X L 1  «P1  ,US1.XN5:  ,Y,XLSDT 

100  rONT[HUc 

if  iv.oRrn)  ioip4*ioi 

101  RETURN 
FNO 

StBFTC  ATP  R94iXR7 

C  SLBROIJTI^  ATP  -  ALTITUDE  T  EMP  A.N[J  PRESS  I  0 1 0 55 / RAMSA V 1 

SL-SROl  Tl^f  ATPUTfST*  jTFST  , KT E ST *L T E $T , A , T4 , PA *  DA . CSA ] 

DIMENSION  HQ! 221 iTRB<22 ) » ELF <22 1 >PBI 221 *FF8<  ?2» 

0 AT  A < HP  111, [  =  1*221/0. *11 300. *2  0000, *32000, ,4  7000, .520  00. *6 1 000 . * 
1790CO. *88743, ,96451 . ,108129. >1 17776. *146541 ,,156071 . , 165571 
21844  65. , 22  1967. ,2 864  76,,  3  76  i  12. ,4  63  526. *&^B2 20, *6  305  30./ 

DATAtTFP I 1 1  i  1  =  1  *22  1  /28  6, 15, 216. 65,2 16. 65. 228, 6 5. 270.65 *270.65. 
1252. 65, 1  80,65  *180. 65. 21 0.65* 260. 65*  360,6$ , 76P.6 5 , ' 110,65.12  10.65, 
21350,65, 1550. 65*1 8 30, 65, 21 £0.65* 24 20, 65, 25 90. 65 *2700.65/ 
f)ATA<cLF<n,I«1.22l/-,O065.C.,,0Cl  ,.<>02  0,0.  *-  ,  flT2  , -.004  .  C  •  * 

1. CC 3C90?3, . 00516636, .U10365 92  . . 02065 66  8  * . 01 57 39 ? 7 , , 0 10526 32 , 

2.00  74  0192  *.0  0 533599, *00  434048, .00  367336, ,00  2981 17 . .00  200699 , 

3,00 ! 33657 ,.00133657/ 

DATAlPPl 1 ) * 1-1, 22 1/101 3.25 *22 5. 32 *64. 74 07, 9, 6301 4 ,1,10905.. 590005* 
1  . 1 02r99  , ,0133777 , 1,64386-03 ,3.0<>75E-0^ .7,3544F-o5,2*5217E-o5, 

25. rfcl7F-^6,3.6943f-C6,2.7926E-06.1 . 6 6 52E-06 ,6 . 9604E-07 , 1 . 5 03 0E- n7 , 
34 , 03 04 r-f  & , l ,69  57E-0fr*3.4  5^2E  C9 *  1 . 1 9 1 e£- 09 / 

HA  TA[FY3(I)*1=L*  22  1/28.9644  ,2  9.9644 , 2  8  .  S644 , 2  0 . 964  <.  ,  7  8 . 9  6  4  4  , 

129, 96 44, 26, 9 64 4, 28. 9644, 26. 564 4, 23, 88, 28,56, 20. 07*26. 92,26. 66* 

226, 4, 25. 85,24.7, 22, *6*39. 94*17. 54, 16. 84, 16. 17/ 

FAC1  =  6.72125 
FAC?  =  0.4700  ? 

CAC*  =  1 545*31 


F AC 5  *  27?1.55S 
A  R  =  P314.32 
G  =  9.80665 

F*0  =  44 

PF  *  6  7  ^  766  ,C 
fFt 1TEST1717, 216*217 
216  ©1=PA*FAC? 

IFlKTCSr-1  )2l3»2f0*2l9 
21B  Pl=Pl»l44. 

219  OC  220  1=2,22 

1FIPBC  n -PI  1222*221 *220 

220  CONTINUE 
721  H* H0 I  1  I 

GO  TO  223 

772  [FI  FIMI  J-U  1235,236*235 

H=  HR  <  !-!  I  -*  <  TMB  I  1-1  l/EtF't  1“1  )  )  *1  iPl/PfM  t~l >  I  *#  «  <  -AR*ELKI  |  )  I  / 

1 <G»FKO) )-l,I 
CO  rrt  22  3 

236  Pf4C  =  P] /PE  r 1-1 ) 

H  =  HB [ I-i  )-JAR»THB4 [-1 ] *A  LOC ( PFAC 1  ] /<G*FHO) 

2=1 PF^Hl/lRE-H) 

225  A*Z/,304fl 
GO  TO  ?r? 

2  17  2=A*,-L048 

?9C  h=:pe*2>/ire-zi 

702  00  ! CO  J=?,Z? 

[FIH-H41  J|  :  101 .1C2.10<- 
1C0  CONTINUE 
1  PI  J=j-1 

1C;  Jl  +  FLMI  Jl  > 

|F< J-?2 )231 * 2^0*230 
??C  ?2J 

GO  T0  2*1  ? 

?■*  1  Ffc'i  =  EHR  I  Jl  -  <  EHB  t  J  T  -FMP  I  J+  1  H  *  l  I  H  HB  I  J  I  )  /  <  H9 f  J*  1  J  -HB I  J  )  )  1 
232  TA  =  <  E*/E”C)*TA 

CSA-=FAC1  *SQRT  <  A  R*  T  A  ) 

FFtFLM  J)  I  103.1^4 ,1C3 

]:«  P  =  5AIJ)M  T4?l  j)/i  IMP  C  J|  *ELM<  J  )*  CH-H1 1  J)  )  »  I*  *  I  <G«MQ  I  /  l  J  >  )  I 

GO  "O  1 f  5 

1 C 4  B  =  PRIJ)*FXP  It - IG*F"U )* I M-Ha t  '  AR-TNBC Jl I ) 

10*=  r,rt  TO  I  106  I  1C71  *  JTEST 
1^6  TA^I ,P»T  A 
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107  1F< lUST I210*21Cp21A 
2|S  6 C  TO  C 10R« 109;*KTE5T 
iCe  pa=P*, £14^03766 
60  TO  2 10 

1T9  FA=P*2. 0865*36 

210  IF  I JTE5T-1  J211»2tl*21? 

711  FAC35-FAC3 
CO  TO  213 
?\2  FAO^fAC^ 

?)3  0A  =  (FM  /FAC3^)*<P*2,O08tS438/TA  | 

CG  To  l  lin, 214*1141 lUTEST 
3 1 A  FAC2 -0  « 

AR=0. 

0=0. 

«-MC=0. 

rf-:. 

F4C1=C, 

FAC  3  =  0* 

FAC 5  =  0  * 

21*.  DA=DA.'32*17404$ 

1 1 0  RCTURN 
FND 

SIB^TC  SIC#/  M94.XR7 

SU^ROUT  INF  SLOW  I XX *2 1 H  ,  J]  , ! T  ,  W,N£RR  3 
C  1  ACF  IS  WRITTEN  WITH  LIVE-  OF  CONST  AM  XX 

C  21111  AND  XX  ARE  INDEPENDENT  VARIABLES 

C  Z  C  J 1  1  is  The  dependent  VARIABLE 

C  AK  =  41,  IF  XX  INCREASES  MONO TON1CALLY  ON  TAPE 

C  AK-  -1,  IF  XX  DECREASES  MCNCTSMCALLY  on  tape 

C  IT*  TAPF  UNIT 

C  NV  =  NO,  OF  VARIABLES  ufv  TAPE  C0R  EACH  XX  .<NOT  GREATER  THAN  9> 

C  NO*  OF  POINTS  FOR  EACH  XX  NOT  GREATER  THA.N  15<" 

C  PFMS  EXECUTION 

DIMENSION  X  U)  »Y  {<.♦<?♦  1  SOI  i  Z  (  9  )  *U<  A  »  •  VI 4  )  *W  C  4  \  *NP  f  4  } 

COMMON  1  T  I  ?U  ) 

I F I IMFTf 11)  7*1 .7 

1  RACK  spacf  n 

RFADI IT)  DUM 
REWIND  IT 

no  2  K«1  *3 

Rr AD f  1  T  )  X<Jt>*J*HY<<*|  pUi  r*l.NV»>L*J*J  ) 


2  NPin  =  - 

XV*X  (  7  )  -X  ( 1  1 
AK  *  AB&<XW>/XW 

om=i, 

I  MF  T  ( 1  )  ■  1 

xx/=xx 

?M  =  3 

GC  TO  70 
7  NFt?R  =  0 

C  EXCEPT  FOR  first  Time  THROUGH 

I f i  (xx-xrwi  ]  >♦  i xx— x i h? > i  noo*ioo»io 
10  T^«p- rxx-xxx 1 *AK 

DTR2-  Af(S(  TEHP  I /TEMP 
GC=  DIR!»D1R2 
XXX -XX 

01*1 =DIR2 
IFID1R?]?0*8*50 

C  NEGATIVE  DIRECTION 

20  lFIC>Ol30ia>40 
30  PACK  SPACE  IT 
PACK  SPACE  IT 
BACK  SPACE  IT 
GO  TO  402 

40  IV  1M-1 

[FIIM  1401.401*402 
401  JK=4 
t*n2  m:  =  jm*i 

BACK  SPACE  IT 
BACK  SPACE  IT 
IF»M1-4>4D4. 404*403 

403  Ml =1 

404  N2=M 14  : 

1 F tf2-4  )406*4Cfi*40S 
403  M2  =»  1 

4  06  REACT  I T  |  XI 1M] ,J,C IY!  IM*I *L>*I=1 *NVt *L=1 *J  ) 
NCH  1M  IrJ 

JFf IXX-X I  Ml) l«(XX-X112|t HOC* 1  DC. 42 
t2  IFtXfMl l-XiM2))A0. 43.40 
C  FRROR.  VARIABLE  OFF  FRONT  END  OF  TAPE 

43  C0NT1NME 
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NERR-l 
GO  TO  200 
C 

C  POSITIVE  DIRECTION 
50  I F  t  GO ) 60  *8 1 70 
60  READ \ IT ]  OUM 
REAOMT)  DUM 
REACH  IT  ]  OUM 
GO  TO  702 
70  1  M=  1  1 

„  IFMH-A>702»702,701 

TOt  1  M=  ] 

702  M 1 = ] M  —  1 

IF<M1)7C3*?03»7C4 
70S  Ml-A 

704  M2  =  M1 -1 

IFIM2I 705*705  *706 

705  M2-4 

706  READIITJ  MtN|fJ»||Y(]M,l»UiW*«V)»L.l,J  I 
NP I f  M ) ■ J 

rFitxx-x<Mi  ij*ikx-k(m?)]|  iDe»iro*?e 

c 

C  TAPE  SEARCH  COMPLETE  *  00  CROSS  FOUR  POINT 

]C0  DO  150  1 *4 

NP*  =  NPUC  1-1 
DO  ns  I «  1  »NPK 

IF<  <Y  (X,  ]  ]  11  Y(K»J1  »t+l  >~Z'  UH  )  12S«1?S»11* 

1  15  CONTINUE 
NE9R-1 
GO  TO  20C 

125  IFt r-! 1127*126il27 

126  J=0 

GO  TO  13? 

127  IF1Z-NPH1 156*135 *136 

135  J=NPK-3 
GO  TO  137 

136  J»J-2 

137  CO  HC  L “  1  .4 
MX=L+J 

UTL  >  =  YIKtl liMXt 
140  VIL  \  s  YIKiJl  »V,X) 


150  CALL  IMTRP  r  4  *  U  » V  »  2  I  I  1  )  » W  I  K  )  ) 

CALL  INI RP  14iX»W»XX»2 (Jl)  1 
175  RETURN 
6  CONTINUE 

200 

2  F  C 1  MET  <  2  I  I  2  I  5  » ?  0 1 »  ?  1  5 
201  ]M£T  <21*1 

00  210  J  M= 1 • 4 
VR1 TE(6*202) XI  INI 
702  FORMAT ( 1 H] Cl 6 >0 > 

?C3  FORma T  f ] h  7E 16  »8 l 
NXXXX=NR{ IM) 

210  WRf TE<6>203H  tYf  >«♦ | *L> »I-1 »NV J  *L-1 tNXXXXl 
215  RETURN 
END 

ilPS’TC  TAUSsS  M9  4  i  XR  7 

SUPROLTTNF  GAUSS  ID*X1 

C  GAUSS  CONTAINS  EIGHT  HLACE iSIXTEEN  POINT  INTEGRATION  CONSTANTS 

DIMENSION  fil 161 .XC16J 
X( ) ) ■«  09SC1 ?5 ) 0 
X C 2  J -.29160355 
XO)». 45601676 

X  C  4  )  =  , £17876  2 4 
X 1 51 =.  7554044  1 
XI 61 =,36563120 
XI 7I=.94457SC2 

XI  6  I =,93940091 
P  I  1  )  =  «  ]  8  945061 
9I?1=.1«26334Z 
0H  1  =  .16915657 

n  I  4  )=, 349*9599 
0151=. 12462897 
BI61=. 09515351? 

3l7l=. 062253524 
0181= .027 15245° 

DO  12145  1=9*16 
J  = 1 7-  I 
XI  I  I =- X  I J  t 
12344  ph  i  =  pi  ji 
°E  T  URM 
END 

^ I P  FTC  INTRD  ^94, XR7 
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_ sum  - 1 . _ _ 

_ _ _  5'JMDbI  ,  _ 

_ DO  201  J«t»N _ 

_ LEiJrIx_20O.2Ql«2CO 

_  2  00  SUMNaSlMN*  I  XINT-Xf  J )  i 

—  '  SlMn«Sl.M!)»IX,  Ll-XUl  ] _ 

_  .  201  CCNTIHIE _ _ 

2P2  YtNT-=vifll4Ym*5J*N'£UfO_ 

_ sxmia _ _ 

CnO 
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